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ABSTRACT 

We show that short-term perturbations among massive planets in multiple 
planet systems can result in radial velocity variations of the central star which 
differ substantially from velocity variations derived assuming the planets are 
executing independent Keplerian motions. We discuss two fitting methods 
which can lead to an improved dynamical description of multiple planet systems. 
In the first method, the osculating orbital elements are determined via a 
Levenberg-Marquardt minimization scheme driving an N-body integrator. The 
second method is an improved analytic model in which orbital elements are 
allowed to vary according to a simple model for resonant interactions between 
the planets. Both of these methods can determine the true masses for the planets 
by eliminating the sini degeneracy inherent in fits that assume independent 
Keplerian motions. We apply our fitting methods to the GJ 876 radial velocity 
data (Marcy et al. 2001), and argue that the mass factors for the two planets 
are likely in the 1.25-2.0 range. 

Subject headings: stars: planetary systems 
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1. Introduction 

About two thousand nearby stars are now being surveyed for periodic radial velocity 
variations which indicate the presence of extrasolar planets (Marcy, Cochran & Mayor 
2000). These search programs have been highly successful, and to date they have discovered 
almost sixty extrasolar planets. 

Recently, systems with more than one planet have been found, and four {v Andromedae, 
GJ 876, HD 83443, and HD 168443) are now known. GJ 876 (Marcy et al 2001) provides 
an especially interesting case. In this system, a combined, two-Keplerian fit to the radial 
velocity data (see Table 1), suggests that the star is accompanied by two planets on orbits 
having a nearly commensurate 2:1 period ratio. The amplitudes of the star's radial velocity 
variations suggest minimum masses of 0.56 Mjup for the inner planet, and 1.89 Mjup for 
the outer planet. GJ 876, an M dwarf star with an estimated mass of O.32±O.O5M (Marcy 
et al 2001), is the lowest mass star known to harbor planets. 

For these orbital parameters, the mutual perturbations of the two planets in the system 
are considerable. In the case of fixed Keplerian ellipses, the orbital elements Pi j2 , ei 2 , T 1)2 , 
u>i t 2, and K 1>2 , are constant, and these, along with the mass of the star, serve to determine 
completely the positions and velocities of the planets at any time in the coplanar case. 
When the interactions between the planets are non-negligible, however, the orbital elements 
change continuously, and so one must also specify an initial epoch in order to determine the 
motion at future times. Every starting epoch corresponds to a different resultant three-body 
motion. 

In Figure 1, we demonstrate this effect by comparing a synthetic radial velocity curve 
generated by the orbital elements given in Marcy et al. (2001) with one generated by an 
N-body integration using the same elements as an initial condition. We used an epoch of 
JD 2450106.2, corresponding to the reported time of perihelion passage for the outer planet. 
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The red line in Figure 1 shows the radial velocity curve which results from the superposition 
of the two Keplerian reflex motions. The black line shows the radial velocity curve resulting 
from the full three-body integration. After three orbits of the outer planet, the motion 
begins to deviate noticeably from the dual-Keplerian approximation. After several years, 
the motions have diverged completely. 

The rest of this paper is organized as follows: in §2 we show that self-consistent radial 
velocity curves are required for systems such as GJ 876. In §3, we derive improved, fully 
self-consistent dynamical fits to the observed radial velocities of the GJ 876 system. In §4, 
we show how dual-Keplerian fits can be improved using an approximate analytic model for 
the interactions between two massive planets in resonant systems. Further applications are 
discussed in §5, which also serves as a conclusion. 

2. Dual Keplerian vs. Self-Consistent Fits 

The orbital elements given in Table 1 (taken from Marcy et al. 2001, based on 54 Keck 
and 16 Lick observations) were derived under the assumption that they are constants of the 
motion. However, for a system such as GJ 876, where the mutual planetary interactions are 
strong, the elements will change quite rapidly on observable timescales. We can therefore 
regard the parameters in Table 1 as a set of instantaneous "osculating" elements, which, 
given a particular starting epoch, correspond to a uniquely determined initial condition. 

Even with the assumption that sin? = 1 for both planets, the orbital parameters in 
Table 1 lead to a wide variety of evolutions depending on the choice of epoch. For some 
starting epochs, the planets are not in the 2:1 resonance, and the system experiences severe 
dynamical instabilities within five years. For other starting epochs, the planets undergo 
librations about the resonance, and the system is stable over long timescales. One can thus 
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ask: are there any starting epochs for which the osculating elements in Table 1 generate an 
evolution which is consistent with the observed reflex velocity of the star? 

We have computed 10,000 synthetic radial velocity curves using N-body integrations 
of the osculating elements for the fit to the Keck data. Each integration used a different 
starting epoch spaced one day apart. The synthetic radial velocities were given a uniform 
offset in order to match the first Keck velocity point obtained by Marcy et al (2001). We 
then compared each synthetic curve to the remaining 53 Keck velocity observations, and 
computed a reduced Xv 2 statistic for the fit. The best fit occurred for an integration with a 
starting epoch of JD 2450722.82. The reduced Xv 2 value for this fit is 12.14, and the rms 
scatter of the velocities about the curve is 64 ms" 1 . Given that the observational errors lie 
in the range 3-5 ms" 1 , this degree of scatter indicates that the best dual-Keplerian fit is a 
poor match to the observed velocities when mutual planetary perturbations are taken into 
account. 

3. A Self-consistent N-body Minimization Scheme 

It is essential to include mutual planetary perturbations when making fits to velocity 
observations of planetary systems resembling GJ 876. One can make a fully self-consistent 
fit which employs N-body integrations to produce a synthetic reflex velocity curve for the 
central star. Starting with the best dual-Keplerian fit to the combined Keck and Lick 
radial velocity data sets, we have used a Levenberg-Marquardt algorithm (Press et al 1992) 
to iterate an improvement to the osculating orbital elements reported in Table 1. Our 
algorithm examines how the x 2 value of the fit depends on variations of all 10 orbital 
elements, and then locates a set of elements for which the reduced Xv 2 fit is at a local 
minimum. We assume that the planets are in co-planar orbits, and we allow the overall 
inclination angle of the system to vary. 
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We have found a self-consistent model for the radial velocity data which has a reduced 
Xu 2 value of 1.46 and an rms scatter of 13.95 ms -1 . This fit is shown in Figure 2 (black 
curve), and represents a significant improvement to the best dual-Keplerian fit reported by 
Marcy et al (2001) which had x^ 2 =l-88. The osculating elements are given in Table 2. We 
note that the combined radial velocity data favors a co-planar inclination of shu = 0.775 for 
the system. There is, however, a broad minimum around this best-fit value. At shu = 0.55, 
the x^ 2 =l-47, and the rms scatter is 14.69 ms" 1 . Only for sini < 0.50 does the fit begin to 
worsen significantly. 

The combined Lick and Keck data sets comprise our best overall view of the GJ876 
system. The 13.95 ms -1 scatter of our best self-consistent fit to this data is due primarily 
to the large uncertainties in the Lick velocities. It is therefore of interest to independantly 
examine self-consistent fits to the Keck data only. Applying our Levenberg-Marquardt 
code to the 54 Keck velocities, we obtain a best fit with an rms scatter of 6.86 ms -1 , and 
X^ 2 =l-59. In this case, the minimum x 2 value occurs for an inclination sini = 0.55, and 
the fit displays a stronger minimum at this particular value. The orbital elements are given 
in Table 2. 

4. Improved Analytic Approximations 

The Levenberg-Marquardt minimization method rapidly locates a dynamical fit to the 
radial velocity data for which the Xv 2 value is at a local minimum. There is no guarantee, 
however, that this minimum is global. This has motivated us to develop an independant 
analytic model for fitting resonant systems. Agreement between the two methods can thus 
increase our confidence that the true configuration of the system has been located. 

We start with a dual-Keplerian model using Jacobi coordinates, where the inner planet 
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moves on an orbit about the star and the outer planet moves on an orbit around the center 
of mass of the inner two bodies. We assume that the planets are undergoing sinusoidal 
librations about the 2:1 mean motion resonance in antiphase to each other. The semi-major 
axes a are thus given by: 



where the amplitude A 1; and initial phase res are treated as model parameters, in addition 
to the mean semi-major axis a 2 of the outer planet. Conservation of energy requires that 



while the resonance frequency, n res is given approximately by n^ es ~ 3ein1(m2/m )a\f c [(a)\ 
for the 2:1 resonance in the case where m 2 is significantly larger than m x (Murray and 
Dermott 1999, Eqs. 8.47 and 8.32). Here a\fa(a)\ ~ 0.750 is a constant that depends only 
on the resonance involved. 

The average mean motions fi\ and n 2 of the planets are related by the condition that 
the rate of change of the resonance critical argument is zero at exact resonance. This 
implies that fix — 2n 2 + — 7r 2 ) = 0, where Ti\ and 7r 2 are the longitudes of periastron 
of each planet. The initial mean anomaly M 2 (0) of the outer planet is treated as a model 
parameter, and the initial mean anomaly of the inner planet Mi(0) is then given by the 
critical argument for the resonance 



ai = ai[l + Ai cos (n rcs t + rcs )] 



a 2 = a 2 [1 - A 2 cos (n res t + res )] 



(1) 




(2) 



a = Mi — 2M 2 + 2(tti - vr 2 ) 



(3) 



where the initial value of a is constrained by the initial values of a\ and a 2 . 



At time t, the mean anomaly of body i is given by 




(4) 
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These integrals are straightforward to evaluate since for each planet a, and hence n, is an 
analytic function of time. 

In this system, the mutual planetary perturbations are sufficiently strong that the 
longitudes of periastron will precess rapidly. We model this by allowing each periastron 
longitude to vary linearly with time, where the precession rates 7Ti and 7r 2 are additional 
free parameters. In practice, we checked that the precession rates derived by the model 
approximately matched those from an N-body integration using the same initial conditions. 
Finally, the orbits of the planets are again assumed to be co-planar, but sini of this plane 
is included as a parameter. 

We used the analytic model to generate synthetic radial velocities for the central star 
and compared these with the Keck and Lick observations given in Marcy et al. (2001). We 
initially generated a randomized population of sets of model parameters and then used a 
genetic algorithm to evolve promising sets towards an improved description of the system. 
At each generation, the genetic algorithm evaluates the degree of fit for each parameter set, 
and cross breeds the best members of the population to produce a new generation. 

Figure 2 (red curve) and Table 3 report an analytic model fit for the combined Keck 
and Lick data set generated by the genetic algorithm. The rms scatter in this case is 13.5 
ms _1 , with reduced xl °f 1-47. This represents a substantial improvement on the best 
dual-Keplerian model in Table 1 (rms scatter of 21.9 ms^ 1 ). The fit is remarkably similar 
to that derived independently using the Levenberg-Marquardt algorithm. The osculating 
elements of a model fit generated by the genetic algorithm for the Keck data alone are also 
given in Table 3. This fit has an rms scatter of 7.0 ms" 1 and a reduced xl of 1.67, which 
is comparable to the fit obtained by the Levenberg-Marquardt N-body technique. As with 
the Levenberg-Marquardt scheme, values of sini ranging from 0.8 to 0.5 yield solutions that 
are similarly good fits to the data. 
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5. Discussion 

The most important benefit of self-consistent dynamical fitting techniques for 
multi-planet systems is the ability to break the sini degeneracy and determine the true 
masses of the extrasolar planets. In a few systems, the true masses can be found when a 
planet transits the parent star, but these cases will be rare and confined largely to planets 
with short periods. The techniques described in §3 and §4 can in principle be applied to 
any system containing more than one planet, given a sufficient baseline of observation. (For 
short baselines, our techniques are most applicable to systems having massive short-period 
planets.) Fischer et al 2000 have shown that roughly half of the planetary systems 
found in the Lick Radial survey show evidence of a second companion. Thus we expect 
that numerous additional multi-planet systems will be forthcoming. In the known cases, 
N-body integrations of the v Andromedae and HD 168443 systems indicate that planetary 
interactions are already producing observable deviations from the multiple Keplerian fit. 

Peale (1994) has used a roughly similar analysis to the one described here to determine 
the true masses and inclinations of the planets orbiting pulsar PSR B1257+12. However, in 
that case, a superposition of Keplerian fits provides a good approximation to the observed 
reflex velocity of the pulsar. This stands in marked contrast to cases such as GJ 876, where 
the planetary interactions are an integral component of the overall motion of the star, 
and an analysis based on small perturbations to these Keplerian fits may not necessarily 
succeed. 

Direct integration shows that the systems represented by the orbital parameters in 
Tables 2 and 3 are stable for at least 10 7 years, despite the presence of strong planetary 
interactions. N-body integrations of the fits to the combined Keck plus Lick data show 
smooth librations about the 2:1 mean- motion resonance and also librations about a secular 
resonance in which the apses precess at the same rate (see Figure 3). These resonances help 
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to maintain stability by preventing close encounters between the planets. 

This long-term stability, combined with the remarkable similarity between the fits in 
Tables 2 and 3, derived using entirely indepedent methods, suggests that these orbital 
parameters are close to those of the real GJ 876 system. The main uncertainty remains 
the precise value of sku, with our coplanar fits suggesting values in the range 0.5 to 0.8. 
The planets may also possess a mutual inclination. We have not modeled this possibility. 
Extending the baseline of observations will make it possible to further refine the orbital 
fit. This, however, may take some time. For example, Figure 4 shows the extent to which 
different fits to the Keck data produced by the Levenberg-Marquardt scheme will diverge 
in future. The black and red curves show the deviation of fits with sini = 0.5 and 0.7 from 
the best fit with sini = 0.55. The differences between these fits will remain less than the 
3cr ~ 15 ms _1 observational errors for several years to come. 
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Figure Captions 



Fig. 1. — Synthetic radial velocity variations for GJ 876 assuming a superposition of 2 fixed 
Keplerian motions for planets with elements given in Table 1 (red line), and an N-body 
integration using the same elements (black line). 

Fig. 2. — Synthetic radial velocity curves for the best fits generated by the Levenberg- 
Marquardt scheme (black) and the improved analytic model (red) for the combined Keck 
plus Lick data. The blue circles show the observations. 

Fig. 3. — The evolution of the critical arguments for the 2:1 mean-motion (solid curve) 
and secular (dotted curve) resonances obtained using N-body integrations of the orbital 
parameters in Table 2 (black) and Table 3 (red) for the combined Keck plus Lick data. 

Fig. 4. — The future divergence of fits to the Keck data by the Levenberg-Marquardt scheme 
for different values of sini The horizontal black line shows the best fit (sin? = 0.55). 
The black and red curves show the deviation from this fit for fits with sini = 0.5 and 0.7 
respectively. The horizontal dotted lines show typical 3a = 15 ms _1 observational errors. 



Parameter 



Inner Outer Inner Outer 
(Keck) (Keck & Lick) 



Period P (day) 


30.11 


61.02 


30.12 


61.02 


K (ms^ 1 ) 


81 


211 


81 


210 


Eccentricity e 


0.29 


0.11 


0.27 


0.10 


u (deg) 


329 


328 


330 


333 


Periastron Time T (JD) 


2450031.4 


2450105.8 


2450031.4 


2450106.2 



Table 1: Best-fit dual-Keplerian elements for GJ876 (from Marcy et al. 2001) 
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Parameter Inner Outer Inner Outer 

(Keck) (Keck & Lick) 



Period (day) 


29.995 


62.092 


30.569 


60.128 


K (ms" 1 ) 


85.256 


203.16 


84.931 


205.351 


Mass (Mj) 


1.06 


3.39 


0.766 


2.403 


a (AU) 


0.1294 


0.2108 


0.1309 


0.2061 


Eccentricity 


0.314 


0.051 


0.244 


0.039 


u (deg) 


51.8 


40.0 


159.1 


163.3 


Mean anomaly (deg) 


289 


340 


356 


173 


Epoch (JD) 


2450602.0931 


2449679.6326 


offset (ms -1 ) 






106.38 


sini 


0.55 


0.78 



Table 2: Osculating orbital elements derived by Levenberg-Marquardt N-body 
integration scheme. 
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Parameter Inner Outer Inner Outer 

(Keck) (Keck & Lick) 



Mass (Mj) 


0.861 


3.16 


0.920 3.08 


a (AU) 


0.1290 


0.2070 


0.1291 0.2067 


Eccentricity 


0.258 


0.054 


0.252 0.046 


oj (cleg) 


206 


217 


210 199 


Mean anomaly (cleg) 


254 


296 


248 314 


Epoch (JD) 


2450000.0 


2450000.0 


sini 


0.60 


0.61 



Table 3: Osculating astrocentric elements derived by the analytic model plus 
genetic algorithm scheme. 
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